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Abstract 

Atmospheric muons play an important role in underwater /ice neutrino 
detectors. In this paper, a parameterisation of the flux of single and mul- 
tiple muon events, their lateral distribution and of their energy spectrum 
is presented. The kinematics parameters were modelled starting from a 
full Monte Carlo simulation of the interaction of primary cosmic rays with 
atmospheric nuclei; secondary muons reaching the sea level were propa- 
gated in the deep water. The parametric formulas are valid for a vertical 
depth of 1.5 5 km w.e. and up to 85° for the zenith angle, and can be 
used as input for a fast simulation of atmospheric muons in underwater /ice 
detectors. 

1 Introduction 



The realization of a kru? scale detector for astrophysical neutrinos is today con- 
sidered one of the most important challenges of the next decade for high energy 
physics and astrophysics [1]. Up to now, reduced scale detectors and prototypes 
have demonstrated the feasibility of the Cherenkov detection technique, while 
larger projects are under way [2]. 

The study of muon bundles in deep underground experiments was originally 
motivated because it could provide some clues on primary cosmic ray flux and 
composition in the PeV region [3]. For neutrino astronomy experiments, where 
neutrino-induced muons are discriminated by selecting upward going tracks, at- 
mospheric muons can represent a major background source. Downward- going 
muons can incorrectly be reconstructed as upward-going particles and mimic high 
energy neutrino interactions; muons in bundles seem to be particularly dangerous. 
Simulations [4] show that a large fraction (at least 40%) of wrongly reconstructed 
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upward going events are induced by muons in bundles. These muons are expected 
to arrive almost at the same time in the plane perpendicular to the shower axis. 
Experimental data [5] show that the arrival time spread of muons in bundles is 
smaller than few ns. 

Direct measurements of the muon flux at sea level exist, as well as paramet- 
ric formulas [6]. On the contrary, atmospheric muon flux in the deep water is 
only roughly known experimentally [7]. Some parameterisations for the energy 
spectrum of underwater single muons versus the vertical depth h and the zenith 
angle are available in literature [8, 9, 10]. These calculations are based on the 
parametric description of the muon flux at the sea level, and use semi-analytic 
methods to solve the muon transport equation in a homogeneous medium. In 
any case, all authors assume a flux of single muons only. 

In this paper, for the first time, parametric formulas are given to calculate 
the flux of muon bundles, taking into account the muon multiplicity and the 
muon energy spectrum in a bundle, as a function of the distance from the shower 
axis. The range of validity extends from 1.5 km down to 5.0 km of water vertical 
depth, and from 0° up to 85° for the zenith angle. A full Monte Carlo simulation 
(section 2) was used to evaluate all the kinematics parameters of muons arriving 
at seven depths. The muon energy spectrum depends on the water/ice vertical 
depth on the zenith angle 6, on the bundle multiplicity m and, for each muon 
in the bundle, on the distance of the muon with respect to the shower axis. The 
paramctcrisation provides a self-consistent generator which can be used as input 
by undcrwatcr/ice experiments for the simulation of the atmospheric muon flux. 
In the equations reported in this paper, the following units are used: depth h 
in kmw.e. — 10^ kg m~^; zenith angle 9 in radians; muon energies E^j^ in TeV; 
radial distances Rin m. 

The plan of the work paper is as follows. In sec. 3 the expression to evaluate 
the flux of muon bundles of any multiplicity is presented. The depth-intensity 
relation for the vertical direction, as well as the dependence of the flux with 
the zenith angle compared with experimental data, is discussed in sec. 4. Since 
the largest fraction of events are single muons, in sec. 5 the differential energy 
spectrum for single muon events as a function of the vertical depth and zenith 
angle is obtained. For muons arriving at the detector in bundles, the energy 
distribution depends also on the muon multiplicity and on the distance of each 
muon from the shower axis. This effect is studied in sec. 6, as well as the 
parameterisation of the muon lateral spread function. In sec. 7, the differential 
energy spectrum of muons in bundles is parameterised in terms of all the four 
variables (muon multiplicity, vertical depth, zenith angle and radial distance). 
The comparison of expected results with the single and double muon average 
energy measured at different depths with the MACRO underground experiment 
is reported in sec. 8. 
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2 The Monte Carlo simulation 



The parameterisation of the multiple muon flux and energy spectra presented in 
this work rehes on a full Monte Carlo simulation of primary Cosmic Rays (CR) 
interactions and shower propagation in the atmosphere (the HEMAS code, sec. 
2.1). The adopted primary CR flux and composition model is described in sec. 
2.2. To optimise the showers production, the primary CRs have been divided 
in flve mass groups, six primary energy intervals and two zenith angle regions. 
Each of these 60 files corresponds to a specific livetime for a 1 km^ detector (the 
livetime spans from few hours for the low-energy bins to few years for the high 
energy ones). The final results have been normalized to the same livetime. 

The muons from the decay of secondary mesons reaching the sea level are 
then propagated down to 5 km of water using the MUSIC code (sec. 2.3). All 
the information concerning muons reaching the seven depths from 2.0 km w.e. 
down to 5.0 km w.e., in steps of 0.5 km w.e., and of their primary CR parents 
are kept in ROOT files. The energy spectrum of muons depends, a part on the 
vertical depth h, on the zenith angle 6, on the muon multiplicity in the shower 
m and on the distance of the muon from the shower axis R. 

2.1 The HEMAS code 

For the primary CR interaction and shower propagation in atmosphere, the 
HEMAS (Hadronic, Electromagnetic and Muonic components in Air Showers) 
code is used [11]. 

HEMAS takes into account all the main physical processes occurring in the 
atmosphere: it computes the first interaction point of primary CR on the basis 
of the input cross sections, propagates electromagnetic and hadronic components 
of the showers considering the actual mean free path of particles, takes into 
account the deflection of charged particles by the geomagnetic field. Hadronic 
interactions in the atmosphere are handled with the hadronic interaction code 
DPMJET [12]. This is a model based on the two component Dual Parton Model 
(the hard and soft components). It also contains a detailed description of nuclear 
interactions. The number of nucleon-nucleon interactions is evaluated from the 
Glauber formalism. The mean free path of CR in the atmosphere is related 
to the inelastic cross sections of primary cosmic rays. For primary protons, we 
have Xp-air{g / cm^) — ^.^'^^^(^fc) and similar relations hold for other primaries or 
secondary produced hadrons. 

In this work the last version of the code (called HEMAS-DPM [13]) has been 
used. It includes the implementation of the Earth curvature, allowing to perform 
correct calculations at large zenith angles. Small angle approximations are used 
in the particle transport along the atmosphere, and the reliability of the code is 
restricted to secondary particles with energies E > 500 GeV. Muons with energies 
E < 500 GeV (at sea level) are not followed through water; in any case, their 
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contribution for depths larger than 1.5 km w.e. is very small, and completely 
negligible at 2.0 km w.e. 

The parameterisations are derived from histograms starting at 2.0 km w.e. 
Nevertheless, the validity range of the results can be extended down to ~ 1.5 
kmw.e., since no particular physical processes are expected to dominate at these 
depths, both for the hadro-production in the atmosphere and for the muon trans- 
port in water/ice. 

The so called prompt muons, the component of the secondary cosmic ray flux 
originateds from the decay of charmed mesons and other short-lived particles 
produced in the interactions of CR with the atmosphere, are not included in this 
simulation. The energy at which the contribution of prompt muons to the sea 
level flux becomes equal to that of muons from tt, K decays is expected to be 
~ 10 TeV to ~ 10^ TeV, depending on the charm production model [14]. 

In the last few years, many other CR simulation codes became available. In 
particular, the CORSIKA code [15] has become a sort of standard in the cosmic 
ray community. However, in this context, HEMAS has been preferred, since it 
was deeply used and cross-checked with MACRO data. In particular, the exper- 
imental muon multiplicity distribution was studied in order to infer information 
on the primary cosmic ray composition [16]. Moreover, some hadronic interaction 
features in HEMAS, first of all the transverse momentum behaviour of muon par- 
ent mesons at TeV energies, were strongly checked by studying the underground 
muon pair distance distribution (the so-called decoherence distribution [17, 18]). 
Many other data-Monte Carlo comparisons, including both shower development 
and hadronic interaction features, can be found in [18]. An extensive study, at 
the level of the ANTARES detector, of the muon fiux obtained with HEMAS and 
CORSIKA simulated showers has been performed [4] and has shown a general 
good agreement. 

2.2 The primary CR flux and composition 

The model adopted for the primary cosmic ray energy spectrum is described 
in [19]. It is a phenomenological model, named the poly-gonato model, that 
uses recent results from direct and indirect measurements of cosmic rays in a 
wide energy range, from 10 GeV to 1 EeV. The direct observations are used 
to extrapolate the energy spectra of each element to high energies. Then, the 
sum of all individual contributions is compared to all-particle spectra from air 
shower measurements. The cut-off for each clement is assumed proportional to 
Z, starting with the proton component at 4.5 PcV. 

A simplified version of the model has been used, neglecting the elements 
heavier than iron and grouping the remaining ones into five mass groups: i) 
protons (A^l); a) helium (A = 4); in) CNO group (A = 14); iv) Mg (A = 24); 
v) Fe and heavy nuclei group (A = 56). 

The minimum primary energy chosen for the simulation is 1 TeV/nucleon. 
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Energy bins (see text) 
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5 xlO^ 




60° - 85° 




8 X 10^ 
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4 X 10=^ 
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6 X 10^ 


5 X 10^ 


1.2 xlO^ 



Table 1: Number of generated interactions for each primary mass group, primary 
energy range and zenith angle interval. 

For each mass group, the production of muons at the sea-level from primary CR 
has been subdivided into six energy ranges: I) 1 TeV - 20 TeV (for He and CNO 
groups the lower hmits are 4 TeV and 14 TeV, respectively); II) 20 - 200 TeV 
(for Mg the lower limit is 30 TeV and for iron 60 TeV); III) 200 - 2 x 10^ TeV ; 
IV) 2 X 10^ - 2 X 10^ TeV; V) 2 x 10^ - 2 x 10^ TeV; VI) 2 x 10^ - 2 x 10^ TeV. 

Finally, the production of muons has been subdivided in two zenith angle 
intervals: 1) 0° < 6' < 60° and 2) 60° < 6' < 85°. The number of generated 
showers for each primary mass group, primary energy range and zenith angle 
interval is reported in tab. 1. The general expression for the flux of each element 
is given by: 

(7c-7z)Ac 



.E 



z , 



(1) 



where the parameters 7^ and Ez are the absolute flux at 1 TeV/nucleus, 
the spectral index and the cut-off energy (knee), respectively and depend on the 
considered nucleus. 7c and ec, which characterize the change in the spectrum 
at the cut-off energy Ez, are assumed identical for all spectra. In Table 2 the 
numerical values for all groups are listed. At the end of the shower simulation, 
HEMAS provides its own standard output files; only events containing at least 
one muon at sea level are kept. 

2.3 Propagation of muons through water 

Muon propagation through water is performed using the MUSIC (MUon Simu- 
lation Code) program [20]. 
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-4.7 
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Table 2: Parameters (from [19]) of the flux of primary CR mass groups used in 
this work. 

MUSIC is a 3D muon propagation code originally written to describe muon prop- 
agation through rock. It uses recent and accurate cross sections of the muon in- 
teractions with matter and it has been modified in order to describe muon energy 
loss through water/ice. It takes into account energy losses due to bremsstrahlung, 
pair production, inelastic scattering and knock-on electron production, consid- 
ered as stochastic processes if the fraction of the energy lost by muon is larger 
than 10~^. The angular deviation of muons is taken into account in the processes 
of multiple scattering, inelastic scattering and pair production. 

3 The multiplicity distribution of muons in bun- 
dles vs. depth and zenith angle 

The muhiphcity distribution of underground muons was experimentally studied 
with large statistics by the Frejus [21] and the MACRO [22] collaborations. The 
expected multiplicity distribution for a given primary mass and energy is known 
to be a negative binomial (NB) distribution. The observed distribution is a 
convolution of NB distributions, which can be described as a power law. Following 
the Frejus paper, the function: 



Mm;h,d) = ^^^^ with v^- ^ (2) 

^ ' ' ^ m'^e*'^) (1 + A-m) ^ ' 



has been used as parametric formula for the flux of bundles with different number 
of muons m at a given depth h and zenith angle Q. Here v\ and A are free 
parameters, depending on h and Q. The phase space has been divided in 7 values 
of vertical depth h (from 2.0 down to 5.0 km w.e. in steps of 0.5 km w.e.) and 
9 values of zenith angle Q (from 0° up to 80° in steps of 10°). Histograms have 
been flUed with all the muons (single or in a bundle) reaching a given vertical 
depth /i, and within = ±1° (±3° for the last bin, due to statistics reasons) 
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centred with respect to each reference zenith angle (i.e., the histograms for the 
value of ^ = 10° have been filled using all muons with zenith angles 9° < < 11°) 
obtained with the full Monte Carlo simulation. The 63 multiplicity distributions 
have been fitted with eq. 2. 

The parameter A 

The Frejus collaboration [21] found, at a depth of ~ 4850 hg/cm'^, a value 
of Apr = 0.66 X 10~^ for shower multiplicities up to 25. A non zero value of 
A would result in an incorrect increase of the flux for muon multiplicities larger 
than m ~ 1/A. All the fits performed are compatible with a value A = 0. In this 
case, ui — 1/ and only two free parameters must be determined in eq. 2. 

The parameter K 

Eq. 2 shows that the parameter K has a direct physical interpretation: for 

m = 1, it is the flux (in units m~^s~^sr~^) of single muons coming from the 6 
direction, at a vertical depth h. As a function of the vertical depth h and of the 
zenith angle 9, it can be described by the equation: 

$(m = 1; h, 9) = K{h, 9) = Ko{h)cos9 ■ e^^^^>""''^ {7x1-"^ s'^ sr'^) (3) 

At a given zenith angle, the flux decreases with depth and two simple expressions 
for Ko{h) and Ki{h) have been found (the values of fltted constants are reported 
in Table 3): 

Ko{h) = Koa ■ /i^"" (4) 
K,(h)^K,a-h + K,b (5) 



The parameter u 

The fraction of multiple muon flux with respect to the single muon flux de- 
pends on the parameter u, which, for a given vertical depth h, is a function of 
sec9: 

iy{h,9)^uo{h)-e'''^^>'"'^ (6) 

For a fixed zenith angle 9, the parameter u increases with increasing vertical 
depth h as: 

i/oih) ^ uoa ■ h'^ + lyob ■ h + voc (7) 

Mh)^'^ia (8) 

As a comparison, the reported value [21] from the Frejus collaboration at an 
average rock depth of 4.85 km w.e. and averaged over all directions, is = 
4.63 ± 0.11. Equation 6 (at the same water/ice depth) produces the following 
values: 1/ = 3.75,4.4 and 13.0 for 9 = 0°,50° and 80°, respectively. 
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Formula 
for the 


Equat. 


Name 


Value 


Equat. 


Name 


Value 


Flux 


4 


Koa 


7.20 X 10-3 


7 




7.71 X 10-2 


(eq. 2) 


4 




-1.927 


7 


I'Ob 


0.524 




5 




-0.581 


7 


l^Oc 


2.068 




5 




0.034 


8 


l^la 


0.030 










8 




0.470 



Table 3: Value of the 9 constants necessary to completely define the flux of bundles of 
any muon multiplicities m, according to eq. 2. The parameter A in eq. 2 is equal to 
zero, as discussed in the text. 

4 The depth-intensity relation 

The vertical direction 

The depth-intensity relation for the vertical direction can be obtained from 
eq. 2. In literature this function is called /(/i, 0) and, according to the notation 
used here, can be written as: 

I{h, 0) = 5^ m • $(m; h, 0) = K{h, 0) • ^^(l/m^''-^)) (9) 

m 

The series at the end of eq. 9 can be numerically calculated; the numerical value 
corresponds to the fraction of muons arriving in bundles with respect to the single 
ones. Fig. 1 shows the I{h, 0) as obtained with this parameterisation, compared 
with [8, 9], where all showers events are assumed as single muons only. The 
line representing the Okada flux is almost indistinguishable from (9). The seven 
points corresponding to the values obtained with the full Monte Carlo simulation 
for ^ = 0° at different depths are also drawn. Fig. 2 shows the vertical flux of 
bundles /(m; h, 0) for different values of the muon multiplicity. The ratio between 
the number of bundles with multiple muons with respect to single muon events 
is ~ 20% at a vertical depth of 2.0 km w.e and becomes ~ 11% at depths larger 
than 4.0 km w.e. 

Zenith dependence 

Fig. 3 shows a comparison of the predicted muon zenith angle distribution 
at two different depths with the measurement of the AMANDA-II under-ice ex- 
periment [23] and of a module of the NESTOR underwater neutrino telescope 
[24]. The upper dot-dashed line corresponds to the calculation presented here, 
at a vertical depth h — 1.64 km w.e. The AMANDA data (marker points) super- 
imposed to the line have been converted to intensities relative to the underwater 
depths, accounting for the lower ice density {pice — 0.917 gcm~^). 
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Full MC simulation 
■ Our fit 



---- Our fit. only single muons 
. - . - Okada calculation 




4.5 5 5.5 
h (km w.e.) 



Figure 1: Muon vertical intensity versus vertical depth. The dashed line represent the 
result of the parameterisation presented in this work, eq. 9, for single muon events. 
The full line represents the total muon flux, which includes the contribution of bundles 
with any muon multiplicity m. Both were calculated with eq. 2 with the constants of 
Table 3. The total flux is in very close agreement with the dashed-dotted line of the 
Okada calculation [8] (the difference between this work and Okada is -1% at 1.5 km 
w.e. and -7% at 5.0 km w.e.). The prediction from [9] is also reported as a dotted line. 
Both [8] and [9] assume a flux of single muons only. The points represent the values 
obtained with the full Monte Carlo simulation. 



One module of the NESTOR neutrino telescope was recently deployed at a 
depth of 3.8 km w.e. In [24] a measurement of the flux of cosmic ray muons is 
reported as a function of the zenith angle 6 according to I{6) = IoCos°'{6), where 
lo = (9.0 ± 0.7) X 10 



-5 



-2—1 

s sr 



and a = 4.7 ± 0.5. In fig. 3, the full line 
represents this calculation for h = 3.8 km w.e., while the dashed lines indicate 
the 1-sigma error band on the fit parameters of the measured NESTOR angular 
distribution. 



5 Energy spectrum of single muons 

The energy-loss processes for TeV muons are theoretically well understood and 
they are usually expressed as the sum of a contribution due to continuous process 
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m=1 (full MC simulation) 
m=2 (full MC simulation) 
m=3 (full MC simulation) 
m=4 (full MC simulation) 
lines = our parametrization 




4.5 5 5.5 
h (km w.e.) 



Figure 2: Flux /(m; h,0) of bundles from the vertical direction {9 = 0°) for different 
values of the muon multiplicity (m = 1,2,3,4) as a function of the vertical depth h. 
The points represent the values obtained with the full Monte Carlo (the errors are 
inside the symbols). The full lines were computed using eq. 2 with the constant of 
Table 3. 



(a) and a term due to catastrophic energy losses (/?): 

-{^) = « + /3i!,. (10) 

Both the quantities a and f3 depend on the medium and on E^^. The expected 
energy distribution of underground muons, assuming a power-law for the primary 
beam energy, at a depth X = hjcosO is [25]: 

-G- i?,e^^(i-^)[E, + e(l - e-^^)]-^ (11) 



d{logioE^ 



where 7 is the spectral index of the primary beam and e = a/(3. It must be 
emphasized that in this context the quantities e, j3 and 7 are considered as simple 
parameters without a specific physical meaning. Eq. 11 provides the shape of 
the energy spectrum, where the constant G represents a normalization factor: 

G(7, (3, e) = 2.30 ■ (7 - 1) ■ e^^^'^ • e^^-^^^'^ • (1 - e-^-^)(^-^) (12) 

The overall flux of single muons is a multiplicative factor, to be calculated with 
eq. 2. The value of the three parameters 7, /3, e depends both on the vertical 
depth h and on the zenith angle 9. 
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Figure 3: Dot-dashed line: predicted total muon flux as a function of the cosine of 
the zenith angle at the depth oi h = 1.64 km w.e. The data points represent the 
AMANDA-II measurement [23]; the error line is inside the marker point. Pull line: 
predicted flux at the depth of /i = 3.8 km w.e. The dashed lines correspond to the 
1-sigma error on the fit parameters of the angular distribution as measured [24] by the 
NESTOR underwater module. 



As before, the phase space has been divided in 7 values of vertical depth 
h (from 2.0 down to 5.0 km w.e.) and 9 values of zenith angle 6 (from 0° up 
to 80°). 63 histograms have been filled with the energy of single muons; the 
histograms have been normalized to unit area, and fitted with eq. 11. In order 
to simplify the parameterisation, the value of the (3 parameter is fixed to /? = 
0.420 {km w.e.)~^ = 4.2 x 10~^ {hg cm~^)^^. The units of e are TeV, and 7 is 
dimensionless. The resulting muon energy is in TeV. 

The parameter 7 

The parameter 7 depends on the vertical depth h {km w.e.) only: its value is 
about 3.8 at 2.0 km w.e., and decreases to about 3.6 at 5.0 km w.e. 7 has been 
parameterised as (see Table 4): 

7 = 7(/i) = 7o ■ /n(/i) + 7i (13) 



The parameter e 

The parameter e depends on both the zenith angle and the vertical depth. At 
a given depth h it shows a linear dependence on sec6, and it can be parameterised 
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— 2.0 km.w.e. 

— 3.0 km.w.e. 
4.0 km.w.e. 

5.0 km.w.e. 

. I .... I I ... I 
-3 -2.5 -2 -1.5 -1 -0.5 0.5 1 

LogioE,(E^inTeV) 

Figure 4: Differential energy spectra of single muons from the vertical direction at 
various depths {h = 2,3,4 and 5 kmw.e.). The marker points (superimposed to the 
h = 2 km w.e. line) correspond to the values obtained with the full Monte Carlo 
simulation. The lines were computed with eq. 11 (constant of Table 4), which gives 
the normalized shape of the distributions. The flux is obtained using a multiplicative 
factor computed at each different depth with eq. 2. 



10' 



10' 



^r\■ 
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as: 

e^e{h,9)^eo{h)-sec9 + ei{h) (14) 
At a given zenith angle, e increases with increasing depth h as: 

eo{h) = eoa-e'"'-'' (15) 
€i{h) ^ eia ■ h + €u (16) 

Comparison of the parametric formula (2) with the full Monte Carlo 

At the end of the fitting procedure, using the seven constants (table 4) it is 
possible to calculate the single muon energy distribution for any vertical depth 
h and zenith angle 9 in the range of validity. The 63 distributions obtained 
from the full Monte Carlo have been compared with the calculated ones. Each 
distribution has a maximum and two flex points. The energy E^"^^ corresponding 
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Formula 
for the 


Equat. 


Name 


Value 


Equat. 


Name 


Value 


Energy 


11 


/3 


0.420 


15 




0.0304 


spectrum 


13 


7o 


-0.232 


15 




0.359 


(eq. 11) 


13 


71 


3.961 


16 


eia 


-0.0077 










16 


ei6 


0.659 



Table 4: The value of the 7 constants necessary to define the (normalized) energy 
spectrum of single muons, eq. 11. 



to the maximum of the distribution is 

(7 



The values of the energy E^"^ and that of the two flex points obtained with the 
parametric formula eq. 11 differ from the corresponding points obtained with the 
Monte Carlo full simulation less than 3%. 

Fig. 4 shows the differential energy spectra of muons coming from the vertical 
direction at various depths. The lines represent the parameterisation obtained 
with eq. 11; note that, as each integrated energy spectrum is normalized to 1, 
eq. 11 gives only the shape of the distribution. A multiplicative factor from 
eq. 2 must be applied to obtain the flux. As an example, the points from the 
full Monte Carlo simulation are superimposed to the curve computed with the 
parametric formula for h — 2 kmw.e. 

Fig. 5 shows the differential energy spectra for single muons at various zenith 
angles and at a depth of 4.5 kmw.e. The lines are obtained from the normalized 
parameterisation of eq. 11 times the absolute flux obtained with eq. 2. 



6 Muons in bundle: the muon lateral spread 

From theoretical and experimental considerations, it results that, in hadron-air 
interactions, particles are produced in clusters; the number of charged hadrons 
follows a negative binomial distribution, whose characteristics depend on the pri- 
mary energy. The transverse momentum pt of the mesons follows in part an 
exponential- law distribution and in part a power- law distribution [13, 17]; most 
of the energy is concentrated in the very forward region (i.e. near the longitu- 
dinal axis). Muons produced in the decay of secondary mesons and reaching a 
given depth of water follow the energy distribution of the parent mesons. As a 
consequence, in a muon bundle, the most energetic muons are expected to arrive 
closer to the axis shower. Therefore, in order to correctly parameterize the energy 
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Figure 5: Differential energy spectra of single muons at various zenitli angles at a 
depth of 4.5 kmw.e. = 4500 hg cm~^. The continuous lines have been computed with 
eq. 11 and constant of Table 4, which gives the normalized shape of the distributions. 
The flux is obtained using the multiplicative factor computed with eq. 2 at each value 
of the zenith angle. 



of muons in a bundle, the muon radial distance R from the shower axis must be 
taken into account. 

The muon lateral distribution in a plane perpendicular to the shower axis can 
be described [11] as: 

dN _ R , . 

dR-^{R + Ro)- 

The average value of the radial distribution {R) depends on the parameters Rq 
and a: (R) — 2i?o/(Q; — 3). Because of the simpler physical interpretation, 
(R), instead of Rq, is used as fit parameter. In equation 18, C represents the 
normalization factor: C = (a — — 2) ■ i?o^^. The parameters a and Rq 
depend on the vertical depth h, the zenith angle 6 and on the muon multiplicity 
m in the bundle. The {h,9,m) phase space has been divided in 189 cells: the 
usual 7 values of vertical depth, 9 values of zenith angle and 3 multiplicities 
m = 2, 3 and > 3. In the following formulas the variable M is defined as: 

M = m, if m < 3 
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M = 4, if m>A 



(19) 



The radial distance R of each muon in a bundle of a given multiplicity M, 
reaching the seven vertical depth h and in a bin of = 1° (3° for the last bin), 
with respect to one of the 9 reference values of zenith angle, was used to fill the 
histogram corresponding to the (h, 9, M) cell. The 189 distributions have been 
fitted using eq. 18. 

The parameter (R) 

The average value (R) of the radial distribution depends mainly on the vertical 
depth (it decreases when h increases). Then, for a given h, (R) decreases with 
increasing of the muon multiplicity. Finally, (i?) docs not depend on the zenith 
angle 6* up to ~ 50°, then it decreases with increasing 9. (R) (in units of m) is 
parameterised as: 

{R) = p{h, 9, M) = po{M) ■ h"' ■ F{9) (20) 

where: 

Po(M) = poa • M + po6 (21) 
^(^) - e(^-^o|./ + i (22) 

The parameter a 

The parameter a increases with the depth h and, at a given depth, it shows 
a smooth decrease with increasing M: 

a = a{h, M) = ao{M) ■ W''^ (23) 

where: 

ao{M) ^ aoa ■ M + aob (24) 

ai{M) ^ aia ■ M + aib (25) 
The value of all constants is reported in Table 5. 

Fig. 6 shows the normalized lateral distribution of double muons for the 
vertical direction and at different values of the vertical depth h. The average 
value of the lateral distribution decreases when h increases because the most 
energetic muons in the bundle arrive closer to the shower axis. 

Fig. 7 shows the normalized lateral distribution of muons with multiplicity 
m = 2, 3 and > 3 (M = 4) from the vertical direction and at the depth of 
3.5 km w.e. The average value of the lateral distribution decreases when M in- 
creases, because showers with large multiplicity were originated by higher energy 
primary CR parents. 
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Formula 
for the 


Equat. 


Name 


Value 


Equat. 


Name 


Value 


Muon 


21 


POa 


-1.786 


24 




-0.448 


lateral 


21 


Pob 


28.26 


24 


O(0b 


4.969 


spread 


18 


pi 


-1.06 


25 




0.0194 


(eq. 18) 


22 




1.3 


25 


aib 


0.276 




22 


i 


10.4 









Tabic 5: The value of the 9 constants necessary to define the (normahzed) distribution 
of radial distances in bundles of muon multiplicity M, zenith angle 9 and vertical depth 
h, eq. 18. 

Note that for the parameterisation of the energy spectrum of single muons, 
the radial distance of the muon from the shower axis is not considered. However, 
the formulas for the muon lateral spread (eq. 18 to eq. 25) are valid also for the 
case M = 1. 

7 Muons in bundle: energy spectra 

The energy spectrum of muons arriving in bundles is described by the same 
function used for the single muons, eq. 11. For muons in bundles, the parameters^ 
7*, /?* and e* depend, apart on the vertical depth h and the zenith angle 6, on 
the muon bundle multiplicity M and on the radial distance R of the muon from 
the shower axis. 

The (/i, 9, M, R) phase space has been divided in 504 cells: 7 values of vertical 
depth; 4 intervals of zenith angle ( 0° - 20°, 20° - 40°, 40° - 60° and 60° - 80°); 
3 values of muon multiplicity (Af = 2, 3, 4); six intervals of radial distance of the 
muons from the shower axis: — 5 m, 5 — 10 m, 10 — 15 m, 15 — 25 m, 25 — 45 m, 
> 45 m. Each [h, 6, M, R) cell correspond to one histogram which has been filled 
with the value of the muon energy from the full Monte Carlo simulation. The 
504 distributions have been fitted with eq. 11. As in the case of single muons, 
7*, f3* and e* have to be considered as fit parameters without any direct physical 
meaning. The value f3* = 0.420 {km w.e.)"^ = 4.2 x 10^^ {hg cm^^)^^ has been 
used. The unit of e* is TeV, 7* is dimensionless. The resulting muon energy is 
in TeV. 

As for single muons, each ^^^^^^^ ^ distribution has a value correspond- 
ing to the maximum (see eq. 17) and two fiex points. The values of E*'"^"-^, as 
a function of 7*, and e*, depend on the vertical depth h, on the zenith angle 9, 

^To avoid confusion with the symbols used for the single muon case, the three parameters 
are indicated with a (*) 
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R{m) 

Figure 6: Lateral distribution (normalized to unity area) for events with two muons in 
the bundle (M = 2). The lines were computed with eq. 18 and the constant of Table 5 
for the vertical direction {6 = 0°) and 4 different depths: ^ = 2,3,4 and h = 5 km w.e. 



on the muon multiplicity M and on the radial distance R of the muon from the 
shower axis. In particular, keeping constant the three remaining variables, the 
value of K*'"*"^: 

- increases when the zenith angle 9 increases; 

- increases when the multiplicity M increases; 

- increases when the vertical depth h increases, reaching a constant value for 

h > 4.5 km w.e.; 

- decreases when the distance R of the muon with respect to the shower axis 
increases. 

The 504 values of the IoqiqE*'"^"'^ obtained with the parametric formula de- 
scribed below, differ from the corresponding values obtained with the full Monte 
Carlo simulation at most by 4%. 

The pcirameter 7* 

As for the case of single muons, 7* does not depend on the zenith angle 6. For 
small values of i? (i? < 10 m) it increases faster than for large values {R > 10 m), 
where it shows a linear dependence on R: 

7* = 7*(i?, h, M) = a{h) ■ R + h{h, M) • (1 - -e'i^^y^) (26) 

The slope a in eq. 26 does not depend on the muon multiplicity M and it 
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R{m) 

Figure 7: Lateral distribution (normalized to unity area) of muons in bundles with 
muon multiplicities M = 2, 3 and M = 4 (m > 3). The lines were computed with eq. 
18 and the constant of Table 5 for ^ = 0° and the vertical depth h = 3.5 km w.e. 



increases with the vertical depth h. Also the intercept b increases with h, but 
with a different rate for different values of M. Parameters a and b can be described 

as 

a — Qo ■ h + ai (27) 

b{h, M) = bo{M) ■ h + bi{M) (28) 
The dependence of bo and bi on M has been parameterised as: 

bo{M) ^boa-M + bob (29) 

bi{M)^bia-M + bu (30) 

The correction factor for the parameter b in eq. 26 has an additional parameter 
which depends on the depth, q — q{h), as: 

Qih) = qo-h + qi (31) 



The parameter e* 

The parameter e* does not depend on the bundle multiplicity M, and it grows 
linearly with the zenith angle 9 (in radians) as: 

e* = e* {R, h, 9) = c{R, h) ■ 9 + d{R, h) (32) 
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Figure 8: Differential energy spectra of muons in bundles with multiplicity M = 2, 
and five different radial distances R from the shower axis: R = 3, 10, 30, 50 and 70 m. 
The lines were computed with eq. 11 and with the constants of Table 6, assuming the 
vertical direction and the depth oi h = 3.5 km w.e. = 3500 hg cm~'^. The marker 
points superimposed for references to the line R = 3 m were obtained with the full 
Monte Carlo simulation and used for the fitting procedure. 



For a given value of the depth h, both the functions c and d in eq. 32 decrease 
with increasing value of the radial distance R. Thus: 

c{R,h) =co{h)-e'''-^ (33) 

d{R,h) = do{h)-R'^'^''^ (34) 

Ci is a constant, while Cq, do, di depend on the vertical depth h as: 

co{h) = coa - h + cob (35) 
do{h) = doa - h + dob (36) 
di{h) = dia - h + dib (37) 

The value of all the constants necessary to define the energy spectrum of muons 
in bundles is reported in tab. 6. 

Fig. 8 shows the normalized energy distribution of double muons (M=2) 
for five values of the radial distance R of the muons from the shower axis. A 
depth of 3.5 km w.e. and the vertical direction are assumed. The maximum of 
the distribution increases when the distance of the muon from the shower axis 
decreases. 



19 



r nTTTin 1 fl 
± yJL lilLlld 

for the 


r (H 11 fl f 


]\] Q rnp 

1 N dliiC 


\/ 1 1 1 P 


r (H 1 1 fl f 

J_L/U Lid Lj . 


]\j Q 1 VI p 
1 N dliiC 


\/ 1 1 1 P 

V dl LLM/ 


Muon 


27 


ao 


0.0033 


35 




-0.069 


in bundles. 


27 




0.0079 


35 




0.488 


Energy 


29 




0.0407 


33 


Cl 


-0.117 


spectrum 


29 




0.0283 








(eq. 11) 


30 


bla 


-0.312 


36 




-0.398 




30 


bib 


6.124 


36 


dob 


3.955 




31 


Qo 


0.0543 


37 


dla 


0.012 




31 


Qi 


-0.365 


37 


dib 


-0.350 



Table 6: The value of the 15 constants necessary to define the (normalized) distribution 
of muon spectrum, eq. 11, for bundles of muons with multiplicity M = 2,3 and 
M = 4 (m > 3), radial distance R from the shower axis, zenith angle 6 and vertical 
depth h. 

8 A comparison of the single and double muon 
energy spectrum with experimental data 

Few measurements of the average energy of underground muons arc available (see 
[26] and references therein). The MACRO collaboration [27] used a transition 
radiation detector to measure the energy of single and double muons reaching the 
detector, crossing different depths of rock. Measurements performed at different 
depths of rock accessible to the detector at different ranges of zenith angles, due 
to the shape of the Gran Sasso mountain coverage, show a higher average energy 
of muons in bundles with M=2 than for single muons. 

To compare the MACRO data with the parameterisation presented here, the 
differences in energy losses of muons in standard rock and water must be taken 
into account. The coefficient for the ionization and excitation (a in eq. 10) 
depends on the ratio (Z/A), which is assumed equal to 0.5 for standard rock and 
is 10/18 for pure water. An explicit formula for the ionization and excitation 
term of the muon energy losses in standard rock and water can be found in sec. 
4.2 of [7]. The value of the ratio (^) (where aR and aw are the ionization- 
excitation coefficients for standard rock and water, respectively) is 0.858 (0.863) 
for a maximum transferable energy of 1 (100) TeV. The coefficient for radiative 
processes (bremsstrahlung, pair production and muon hadroproduction, [3 in eq. 
10), is larger in standard rock than in water. In the 1-100 TcV muon energy range, 
(3w = (3.15H-3.8) X 10-^ {hg-^cm^) and /5r = (4.0H-4.7) x 10"^ (hg-^cm^), where 
Pr and Pw are the averaged values of the coefficient P in standard rock [25] and 
water [8], respectively. 
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As a consequence, the residual energy of muons crossing a given amount of 
water (in hg/crn^) is higher than the residual energy of muons crossing the same 
amount of standard rock. The value of the average energy -E^ ^ocA: of single and 
double muons crossing a given depth h of rock is estimated^ from the corre- 
sponding value of the average energy E^y^r of single and double muons crossing 
the same depth of water, using the correction factor: 

E,,noc, = (^) • [J _ ■ E,,w = r^(h) . E,,^ (38) 

In Table 7 the values of the average energies reported by the MACRO collab- 
oration at four different depths arc compared with the values obtained with the 
present parameterisation, after the correction for the different media, eq. 38. For 
the comparison, a second effect must be taken into account. Due to the shape of 
the Gran Sasso mountain, a given range of rock depth in MACRO corresponds 
to a range of zenith angles. In fact, for increasing zenith angles, also the average 
rock depth crossed by the muon increases. Using the Gran Sasso rock map [28] 
(rock thickness as a function of zenith and azimuth angles), the average value 
of the zenith angles for muons in each rock depth bin (column 2) is estimated, 
taking into account that only muons with 9 < 45° were accepted for the MACRO 
measurement. For double muons only, the weighted value of the average radial 
distance of the muons from the shower axis, reported in column 3 has been used. 
The computed values of the correction factor T^{h) are reported in column 4; 
the central value of the ratios computed at 1 and 100 TeV is considered. 

The excellent agreement between data and results from the parameterisation 
presented here (see table 7) can be partially biased by the fact that some HEMAS 
features were tuned using MACRO data (although data reported in [27] were not 
considered) . 

9 Discussion and conclusions 

The knowledge of the multiple muon flux is an important requirement for neutrino 
telescopes. The measurement of the muon flux in a homogeneous medium, such 
as sea water or ice, is also of fundamental importance for the study of the muon 
energy losses and to test the flux of primary CR and their interactions in the 
upper atmosphere using the depth-intensity relation. Data obtained in deep 
mines or underground laboratories cannot be used to operate a fine tuning of the 
theoretical formulas, because of the uncertainties in rock composition and density 
along the muon trajectories, as pointed out in [7]. Finally, the measurement of 
the muon flux is important to constraint the flux of high energy atmospheric 

^The formula for the average value for the energy distribution of eq. 11 has the form of eq. 
17, with the replacement of (7 — 1) with (7 — 2) at the denominator. 
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Single muons 


Double muons 


h 

{hg cm~'^) 


e 

(deg) 


R 

(m) 


^ R 


(GeV) 


i R ■ ^lfl,W 

(GeV) 


{GeV) 


J- R ■ J^2,i,W 

{GeV) 


3280 


3 


5.2 


0.768 ±0.02 


250 ± 17 


248 


321 ± 23 


329 


3420 


10 


4.7 


0.765 ± 0.02 


262 ± 18 


254 


366 ± 24 


350 


3600 


20 


4.3 


0.762 ±0.02 


278 ± 19 


265 


400 ± 25 


383 


3800 


30 


4.0 


0.758 ±0.02 


283 ± 19 


278 


417 ±25 


412 



Table 7: Average energy of single (column 5) and double (column 7) muons measured 
by the MACRO collaboration [27] for different values of the slant depth (column 1) and 
of the average muon zenith angle direction (column 2). The predictions of this work 
are reported in column 6 and 8 for m = 1 and m = 2, respectively. The calculated 
values were corrected to take into account the muon energy losses differences in water 
and rock, computed with eq. 38, and reported in column 4. For double muons, the 
computed average distance of muons from the shower axis, reported in column 3, has 
been assumed . 

neutrinos. The absolute normalization of the atmospheric neutrino flux is still 
affected by ~ 30% error in the 100 GeV-1 TeV range. For higher energies, the 
uncertainty can be even larger [29]. 

To give an estimate of the overall uncertainty on the absolute muon flux, 
formulas presented in this work are compared with experimental data and with 
other calculations. An overview of the experimental measurements of the muon 
vertical intensities in the deep sea is in [9], where the data are also compared 
with a parameterisation described there (the comparison with the AMANDA-II 
data is in [23]). The same data set was compared in [7] with the Miyake formula 
[30] , whose predictions are close to the parameterisation of Okada [8] (shown in 
Fig. 1, together with the parameterisation discussed in the present work and the 
Bugaev one). The results presented here differ at most by 7% from the Okada 
parameterisation, and by 15% from the Bugaev one. The Okada and Bugaev 
predictions differ at most by 15%. 

The Monte Carlo code used for simulation (HEMAS) was optimized to repro- 
duce the multiplicity distribution of muons in bundles and the lateral distribution 
of muons inside the bundle at different depth of standard rock. The use of para- 
metric formulas well reproduce all the Monte Carlo distributions. The average 
muon multiplicity and the average value of the radial distribution obtained with 
the full Monte Carlo differ from the results from parametric formulas less than 
5%. The comparison between the average muon multiplicity at different depths 
using HEMAS and CORSIKA agrees at the level of 10%. 

The uncertainty on the energy spectrum and on the average energy of sin- 
gle and multiple muons is more difficult to evaluate. The main uncertainties 
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still arise from the hadronic interaction model and from the muon propagation 
through water. In the distributions resulting from this paramctcrisation the po- 
sition of some fiducial points (maximum of the distribution, flex points) differ 
at most by 4% from the results obtained with the full Monte Carlo. These nu- 
merical approximation can be considered as negligible with respect to the overall 
theoretical uncertainties affecting the full Monte Carlo simulation. To test the 
agreement between the predictions obtained with these parametric formulas and 
experimental data, the average energy of single and double muons measured by 
the MACRO collaboration at different depths [27] have been used. A correction 
procedure has been used to take into account the different media (water in this 
work, standard rock for the MACRO data). All the values agree at level of 5%. 

Finally, it should be noted that the contribution from the so called prompt 
muons is neglected in HEMAS simulation. An unknown uncertainty factor due to 
this process should be included for muon residual energies higher than ~ 10 TeV. 

In conclusion, parametric formulas for the flux of single and multiple muons 
for the interval 1.5 < h < 5.0 km w.e. and for 9 < 85° are given. The energy 
spectrum of single muons and multiple muons in bundles was also parameterised, 
taking into account the dependence of the muon energy on the shower multiplicity 
and on the distance of the muon from the shower axis. The main results can be 
summarized as follows: i) the distribution of the muon multiplicities in a bundle 
does depend on the vertical depth h and zenith angle 9; ii) for a flxed zenith angle 
9, bundles with high multiplicity are suppressed when h increases; in) the average 
distance of each muon in a bundle does depend on the bundle multiplicity m, and 
on h and 9; iv)ai a given depth, the average distance of each muon in a bundle 
decreases slightly when the multiplicity m increases. This can be qualitatively 
understood because bundles with high multiplicities are produced on average 
from primary parents of higher energies; f )the average energy of a muon in a 
bundle depends on the depth h, on the zenith angle 9, on the bundle multiplicity 
m and on the distance R of the muon from the shower axis. 

The parametric formulas presented in this paper can be easily implemented 
in a Monte Carlo generator which can be used to study the response of underwa- 
ter/ice detectors to the flux of atmospheric muons. 

Acknowledgements This work was motivated through studies performed 
with the ANTARES and NEMO Collaborations. In particular, we would grate- 
fully acknowledge the help and the discussions with J. Brunner, G. Carminati, 
S. Cecchini, C. Di Stefano, G. Giacomelli, E. Korolkova and V. Kudryavtsev. 
We also gratefully thank Dr. Paolo Desiati of the University of Wisconsin for 
providing us the reported AMANDA-II data. 



23 



References 

[1] J. Learned and K. Mannheim, Ann. Rev. Nucl. and Part. Science 50 (2000) 
679. 

[2] http://antares.in2p3.fr 
http://www. nestor. org.gr 
http:/ / amanda. uci. edu 
http://www.ifh. de/baikal/baikalhome.html 
http://nemoweb.lns.infn.it/publication.htm 
http:/ / www. kmSnet. org 

[3] T.K. Gaisser, T. Stanev, NIM A235 (1985) 183. 

[4] S.Cecchini at al. (for the ANTARES CoU.), Atmospheric muon background 
in the ANTARES detector, 29th ICRC, Pune, India 

[5] S.P. Ahlen et. al. (The MACRO Collaboration), Nucl. Phys. B370 (1992) 
432. 

[6] S. Miyake, Proc. 13th Inter. Conf. Cosmic Rays, Denver, 5 (1973) 3638. 

[7] P.K.F. Grieder, Cosmic Rays at Earth, Elsevier, Amsterdam, 2001. 

[8] A. Okada, Astrop. Phys. 2 (1994) 393. 

[9] E.V. Bugaev et al., Phys. Rev. D58 (1998) 054001. 

[10] S.I. Klimushin et al., Phys. Rev. D64 (2001) 014016. 

[11] C. Forti et a/., Phys. Rev. D42 (1990) 3668. 

[12] J. Ranft, Phys. Rev. D51 (1995) 64. 

[13] G. Battistoni et al, Astropart. Phys. 3 (1995) 157; INFN/AE-99/07 (1999). 

[14] T.S. Sinegevskays, S.I. Sinegovsky Phys. Rev. D63 (2001) 096004. 

[15] J.N. Capdevielle et al.. the Karlsruhe extensive air shower simulation code 
CORSIKA, KFK Report (1992) 4998; FZKA (1998) 6019. 

[16] MACRO Collaboration, M. Ambrosio et al, Phys. Rev. D56 (1997) 1407; 
Phys. Rev. D56 (1997) 1418. 

[17] MACRO Collaboration, M.Ambrosio et al, Phys. Rev. D60 (1999) 032001. 

[18] M. Sioli, A new approach to the study of high energy muon bundles with the 
MACRO detector at Gran Sasso, Ph.D Thesis, Bologna University, 2000, 
hep-ex/0209029 



24 



[19] J. Horandel, Astrop. Phys. 19 (2003) 193. 

[20] P. Antonioli et al., A three dimensional code for muon propagation through 
the rock, Astrop. Phys. 7(1997)357; arXiv: hep-ph/9 705408. 

[21] C. Berger et. ah, Phys. Rev. D40 (1989) 2163. 

[22] M. Ambrosio et. al. (The MACRO Collaboration), Phys. Rev. D56 (1997) 
1407. 

[23] P. Desiati, K. Bland et al. (AMANDA Collaboration), Response of 
AMANDA-II to Cosmic Ray Muons, 28th ICRC, Tsukuba, Japan, HE 2.3, 
1373-1376 

[24] G. Aggouras et al. (The NESTOR Coll.), Astrop. Phys. 23 (2005) 377. 
[25] P. Lipari, T. Stanev, Phys. Rev. D44 (1991) 3543. 

[26] M. Ambrosio et. al. (The MACRO Collaboration), Astrop. Phys. 10 (1999) 
11. 

[27] M. Ambrosio et. al. (The MACRO Collaboration), Astrop. Phys. 19 (2003) 
313. 

[28] M. Ambrosio et. al. (The MACRO Collaboration), Phys. Rev. D52 (1995) 
3793. 

[29] M. Ambrosio et. al. (The MACRO Collaboration), Eur. Phys. J. C36 (2004) 
323. 

[30] S. Miyake, J. Phys. Soc. Japan, 18(1963) 1093. 



25 



